Thermodynamic stability of Fe/O 
solid solution at inner-core conditions 
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We present a new technique which allows the fully ab initio calculation of the chemical potential of 
a substitutional impurity in a high-temperature crystal, including harmonic and anharmonic lattice 
vibrations. The technique uses the combination of thermodynamic integration and reference models 
developed recently for the ab initio calculation of the free energy of liquids and anharmonic solids. 
We apply the technique to the case of the substitutional oxygen impurity in h.c.p. iron under Earth's 
core conditions, which earlier static ab initio calculations indicated to be thermodynamically very 
unstable. Our results show that entropic effects arising from the large vibrational amplitude of the 
oxygen impurity give a major reduction of the oxygen chemical potential, so that oxygen dissolved 
in h.c.p. iron may be stabilised at concentrations up a few mol % under core conditions. 



The thermodynamic stability of oxygen dissolved in 
iron is a key factor in considering the physics and chem- 
istry of the Earth's core. We present here a new technique 
which allows the ab initio calculation of the chemical po- 
tential of an impurity in a high-temperature solid solu- 
tion, including harmonic and anharmonic lattice vibra- 
tions. We report the application of the technique to sub- 
stitutional oxygen dissolved in hexagonal close-packed 
(h.c.p.) iron at Earth's core conditions, and we show 
that the Fe/O solid solution is thermodynamically far 
more stable than expected from earlier work. The new 
technique should hnd wide application to a range of other 
earth-science problems. 

It has long been recognised that the liquid outer core 
must contain a substantial fraction of light impurities, 
since its density is 6 — 10 % less than that estimated 
for pure liquid Fe 0^]; similar arguments suggest that 
the inner core contains a smaller, but still appreciable 
impurity fraction [gj. The leading impurity candidates 
are S, Si and O, and arguments have been presented for 
and against each of them g]. Ringwood Q and oth- 
ers H have argued strongly on grounds of geochemistry 
that oxygen must account for a large part of the im- 
purity content. However, it has proved difficult to as- 
sess these ideas, because the Fe/O phase diagram is so 
poorly known at Earth's core conditions. (For reference, 
we note that the pressures at the core-mantle bound- 
ary, the inner-core boundary (ICB) and the centre of the 
Earth are 136, 330 and 364 GPa respectively; the tem- 
peratures at the core-mantle boundary and the ICB are 
poorly established, but are believed to be in the region 
of 4000 and 6000 K respectively.) 

The thermodynamic stability of dissolved oxygen is 
governed by the free energy change in the reaction 



(n - l)Fe(solid) + FeO(solid) 



Fe„0 (solid solution) 
(1) 



Let AG be the increase of Gibbs free energy as this 
reaction goes from left to right, excluding the config- 
urational contribution associated with the randomness 
of the lattice sites occupied by dissolved O. Then the 
maximum concentration (number of O atoms per crystal 
lattice site) at which dissolved O is thermodynamically 
stable with respect to precipitation of FeO is c max = 
exp(— AG/fee? 1 ). Several years ago, Sherman |(| used 
ab initio calculations based on density functional theory 
(DFT) to calculate the zero-temperature limit of AG, i.e. 
the enthalpy AH of reaction (1). He found that AH is 
very large (~ 5 eV at the ICB pressure of 330 GPa), and 
concluded that the concentration of dissolved O in the 
inner core must be completely negligible. His argument 
has been widely cited. However, these were static, zero- 
temperature calculations, which entirely ignored entropic 
effects. We shall show here that the high-temperature en- 
tropy of dissolved O produces such a large reduction of 
free energy that Sherman's argument should be treated 
with caution when considering core temperatures. 

Our ab initio calculations are based on the well es- 
tablished DFT methods used in virtually all ab initio 
investigations of solid and liquid Fe including 
Sherman's We employ the generalised gradient ap- 
proximation for exchange-correlation energy, as formu- 
lated by Perdew et al. fllq] , which is known to give 
very accurate results for the low-pressure elastic, vibra- 
tional and magnetic properties of body-centred cubic 
(b.c.c.) Fe, the b.c.c. — * h.c.p. transition pressure, and 
the pressure-volume relation for h.c.p. Fe up to over 
300 GPa @£l). We use the ultra-soft pseudopotential 
implementation |lj| of DFT with plane-wave basis sets, 
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an approach which has been demonstrated to give re- 
sults for solid Fe that are virtually identical to those of 
all-electron DFT methods [[ll| . Our calculations are per- 
formed using the VASP code [B , which is exceptionally 
stable and efficient for metals. The technical details of 
pseudopotentials, plane-wave cut-offs, etc. are the same 
as in our previous work on the Fe/O system p5[. 
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FIG. 1. Gibbs free energy AG of the reaction 
(n — l)Fe(solid) + FeO(solid) — » Fe n O (solid) from ab ini- 
tio calculations. Upper curves show zero-temperature values 
(i.e. enthalpies) for 25 mol % concentration from Refs. Jl^] 
(solid curve) and ^| (long dashes) and for 1.5 % concentra- 
tion (short dashes) from present work. Short lower curves 
show high-temperature results for AG from present work at 
5000 K (dots) and 7000 K (chain line). 

We first report static zero-temperature results for the 
enthalpy AH of the reaction (1). Sherman's results 
later confirmed by the present authors [ ft5"[ , were obtained 
for the high O concentration of 25 mol %, corresponding 
to n = 3, but here we wish to focus on the dilute limit, 
and we take n = 63, which gives a mole fraction of 1.5 %. 
To do this, we treat a 64-atom supercell with the h.c.p. 
structure containing a single O substitutional, and we 
calculate the total ground-state energy and pressure for 
a sequence of atomic volumes, with all atoms relaxed to 
their equilibrium positions at each volume. The enthalpy 
of the pure iron system is obtained from total-energy and 
pressure calculations for a single unit cell of the h.c.p. 
crystal. For the FeO crystal, we obtain the enthalpy from 
total-energy and pressure calculations on a unit cell of 
the NiAs structure. (The high-pressure stable structure 



of FeO is believed to be either NiAs or inverse-NiAs; the 
relative stability of the two structures has been contro- 
versial Jl6|-|l8|,but our own ab initio calculations indicate 
that the NiAs structure is slightly more stable at pres- 
sures above ca. 145 GPa.) The enthalpy AH is reported 
as a function of pressure in Fig. 1, where we also show 
Sherman's results and our own for the 25 mol % case. 
We see that AH at 1.5 mol % is between 0.5 and 1.0 eV 
lower than at 25 mol %, having a value of ca. 4.7 eV at 
330 GPa, but this is still very large and Sherman's ar- 
guments would remain valid if this represented a good 
estimate of the free energy of reaction (1) at Earth's core 
temperatures. 
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FIG. 2. Calculated energy (eV units) as function of dis- 
placement (A units) of O impurity from its equilibrium site 
in basal plane (solid curve) and along c-axis (dashed curve). 
Results are for crystal volume of 6.97 A 3 /atom. 

To get an idea of the freedom of movement of the sub- 
stitutional O atom, and hence its vibrational entropy, 
we now perform a series of calculations in which the O 
atom is displaced by different amounts from its equilib- 
rium site, with all other atoms held fixed at the equi- 
librium positions they have when the O atom is at its 
own equilibrium position. Results for the energy varia- 
tion with displacement along the c-axis and in the basal 
plane for the crystal volume of 6.97 A 3 /atom are shown 
in Fig. 2. Since we are interested in core temperatures, 
we also mark the energy k-gT for T = 6000 K. We see that 
in the energy range set by this T the energy surface is 
extremely anharmonic, with almost vanishing curvature 
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at the equilibrium site and large curvatures for large dis- 
placements. Clearly, the vibrational root-mean-square 
(r.m.s.) O displacement will be ~ 0.45 A (our direct 
ab initio molecular-dynamics calculations confirm this 
value). For comparison, we estimate the r.m.s. O dis- 
placement in FeO at the same P and T to be ~ 0.23 A. 
This implies a substantial vibrational entropy for the O 
impurity, because it fits so loosely into the Fe lattice. 
For the same reason, the vibrations of the 12 Fe neigh- 
bours of the O impurity will be softened, and this will 
also increase the entropy. 

To make quantitative statements about high- 
temperature behaviour, we need to calculate the ab ini- 
tio Gibbs free energies, rather than zero-temperature en- 
thalpies. The ab initio calculation of the Helmholtz free 
energies F of the Fe and FeO perfect crystals is straight- 
forward in the harmonic approximation, since this re- 
quires only the energy of the static lattice and the ab 
initio lattice vibrational frequencies, which we calculate 
by the small-displacement method, as discussed in sev- 
eral previous papers jl0| , |l9| -^H . From F, we then directly 
obtain the Gibbs free energy G = F + PV, by calculating 
the pressure P as —(dF/dV)T, with V the volume. The 
only difficult part of the present problem is therefore the 
calculation of F (and hence G) for the Fe„0 crystal con- 
taining the substitutional O atom. This free energy must 
include the vibrations of many shells of neighbours of the 
O impurity. The harmonic approximation will clearly not 
suffice. We meet this challenge by drawing on recently 
developed ab initio methods for calculating the free en- 
ergies of liquids and anharmonic solids |l(],|22|,^3| . These 
methods rely on two things: empirical reference models, 
parameterised to accurately mimic the ab initio energies; 
and the technique of 'thermodynamic integration', used 
to determine free energy differences. Our overall strategy 
will be to obtain the ab initio free energy F^Lq of the O- 
substitutional system by starting from the ab initio free 
energy F^} of pure Fe and using thermodynamic integra- 
tion to compute the free energy change F^ Q — F^} that 
results from converting a single Fe atom into an O atom. 

We recall briefly that thermodynamic integration 124] 
is a general technique for calculating the difference of 
free energies F\ — Fq of two systems containing the same 
number N of atoms but having different total-energy 
functions Uq(ti, r%, . . , rjy) and Ui(r 1: r 2 , . . . rjv), with 
Yi{i — 1,2... TV) the atomic positions. The technique 
relies on the equivalence between the free energy differ- 
ence and the reversible work done on switching the total 
energy function continuously from Uq to U\. The work 
done is: 

Fx-F Q = f d\(U 1 -U )x , (2) 
Jo 

where the thermal average ( • }\ is evaluated in the canon- 
ical ensemble generated by the switched energy function 
U\ defined by: 



U x = (1 - X)U + \U X . (3) 

To apply this in practice, we use molecular dynamics 
simulation to evaluate the average {U\ — Uq)\ at a se- 
quence of A values and we perform the integration over 
A numerically. 

In principle, we could calculate F^, — Fp^} by iden- 
tifying Uq and U± as the ab initio total energy functions 
and Up*, Q of the pure-Fe and O-substitutional sys- 
tems, but this brute-force approach is computationally 
prohibitive at present. It is also unnecessary, since ex- 
actly the same result can be achieved much more cheaply 
by using empirical reference models. In our recent work 
on liquid Fe [ pd| , we have found that a simple inverse- 
power pair potential <p(r) = A/r a reproduces the ab 
initio total energy very accurately; for the anharmonic 
high-temperature Fe crystal, a linear combination of this 
pair-potential model with an ab initio harmonic descrip- 
tion has been very effective Jl(| . We denote the total en- 
ergy of this latter anharmonic model by C/p° (ri, . . . rjv), 
where rj are the atomic positions. 

In order to make a reference system for the Fe/O sys- 
tem containing a single substitutional O atom, whose to- 
tal energy is Up° l ^ (ri, . . . rjy) (ri is the position of the 

O atom), we simply delete all terms in C/p° f involving 
atom 1 and replace them with a pair interaction poten- 
tial x( r ) = B/r@ . All parts of U F ^ involving the N — 1 
Fe atoms 2, 3, ... N remain precisely as they are in U F C J . 
Our procedure for determining B and (3 starts by requir- 
ing that the pair potential x( r ) should reproduce as well 
as possible the dependence of the total energy on position 
of the O atom with all Fe atoms held fixed, i.e. the curves 
shown in Fig. 2. The values for B and (3 thus obtained 
give us an initial form for U F ^, Q . This initial U F ^, is 
then used in a classical MD simulation to generate a long 
trajectory at the temperature of interest, from which we 
take 100 statistically independent configurations. The 
full ab initio energies are calculated for these configura- 
tions, and the B and parameters are readjusted to give 
a least squares fit to these energies. Finally, the U F °^ Q 
obtained from the new B, (3 is used to generate a further 
100 statistically independent configurations, and B and 
(3 are adjusted once more to fit the ab initio energies of 
these configurations. The B and (3 produced by this final 
step are found to be essentially identical to those in the 
previous step, and we accept them as the optimal values 
for the assumed form of x( r )- 

The free energy difference i*jw — ^Fc 1 between the 
Fe/O and pure Fe systems is now expressed as the sum 
of three contributions: the difference F^, — F0 be- 
tween the reference systems, and the two differences 
Fp^ Q — F F f, Q and F F f - F^ 1 between the ab initio and 
reference systems. All three of these differences are cal- 
culated by thermodynamic integration. We emphasise 
that, although the reference systems play a vital role, 
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the final result for F^, — does not depend on how 
they are chosen. 

Thermodynamic integration to get fp°/ — Fp* { is 
easy and rapid, since only simple potential models are 
involved. We have tested size effects by using simu- 
lated systems containing up to 768 atoms, but we find 
that with only 288 atoms the size-effect errors are less 
than the statistical errors of ca. 30 meV. For the ab- 
initio /reference differences F^, o —F0, o and F^ — F^f, 
the fluctuations of U^, Q — and — Upf are so 

small that explicit thermodynamic integration over A is 
unnecessary, and we can use instead the small- A approx- 
imation explained elsewhere ^5|. We have studied the 
size errors for these ab initio/ 'reference differences, and 
we find that results obtained with 36, 64 and 96 atoms 
are identical within statistical errors. The overall sta- 
tistical error on the ab initio difference F^, Q — F^} is 
ca. 90 meV. We have repeated all the above calculations 
at the four volumes 6.86, 6.97, 7.20 and 7.40 A 3 /atom, 
and from the dependence on volume we obtain the pres- 
sure change on replacing Fe by O and hence the Gibbs 
free energy difference GpL Q — Gpl- 

Our calculated Gibbs free energies AG for reaction (1) 
are displayed in Fig. 1 for the two temperatures 5000 and 
7000 K. We note the very large entropic lowering of AG, 
which, at P = 330 GPa comes down from 4.7 to ca. 1.7 eV 
at the temperature T ~ 6000 K expected at the ICB. 
This is still a substantial positive value, but implies that 
the stability-limit concentration c max = exp(— AG/k^T) 
is ca. 3 mol %, which is far from negligible. 

In assessing our c max value, one should note the re- 
maining uncertainties in our calculations. First, we have 
ignored anharmonicity in the pure Fe and FeO crystals. 
Our recent work on the effect of anharmonicity in pure 
Fe ]2j| showed that at the melting point, anharmonicity 
can contribute as much as 70 meV/atom to the free en- 
ergy; the same might be true of FeO. These effects could 
shift AG by perhaps 0.15 eV. Second, there is the ques- 
tion of strong electronic correlation in FeO, which is a 
prototypical Mott insulator at low pressures. Such corre- 
lation effects will be much weakened at Earth's core pres- 
sures, but could still shift AG by a few tenths of an eV. 
This means that our prediction for c max at a given tem- 
perature is probably not reliable to better than a factor 
of 3. We are therefore cautious about the detailed nu- 
merical value of c max , and claim only that it could be a 
few mol % at ICB pressure and temperature. 

In summary, we conclude that, because substitutional 
oxygen atom fits so loosely into the Fe lattice and has so 
much freedom of movement, it undergoes a very large en- 
tropic lowering of free energy at high temperatures, this 
lowering being as much as 3 eV at 6000 K and 330 GPa. 
The consequence is that substitutional O dissolved in 
h.c.p. Fe may be thermodynamically stabilised at con- 
centrations up to a few mol %. Earlier ab initio calcula- 



tions H which ignored entropic effects should therefore 
not be taken at face value. Finally, we point out that a 
wide range of geological problems depend on an under- 
standing of chemical potentials - for example, all prob- 
lems concerned with the partitioning of elements between 
coexisting phases. The ab initio techniques for calculat- 
ing chemical potentials outlined here should therefore be 
of wide interest. 
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